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1 Introduction 

Blind Iterative Deconvolution (BID) is a technique which was originally developed to 
correct the degrading effects of atmospheric turbulence on astronomical images from single 
short exposure, high signal-to- noise-ratio frames. At the Center for Astrophysics, we have 
implemented a version of BID following the general approach of Ayers and Dainty (1988), but 
extending the technique to use Wiener filtering (Nisenson el al, 1990) and developed it for 
application to high energy images from Einstein and ROSAT. In the optical, the point spread 
function (psf) that degrades the images is due to a combination of telescope and atmospheric 
aberrations. At high energies, the degrading function is the instrument response function, 
which is known to be time and energy level unstable. In both cases the psf is poorly known, 
so BID can be used to extract the psf from the image and then deconvolve the blurred image 
to produce a sharpened image. 

BID follows the general approach of constrained iterative techniques such as those 
developed for phase retrieval (Gerchberg and Saxton, 1972; Fienup, 1978) with blind de- 
convolution (Lane and Bates, 1987). One starts with an image which is degraded by some 
blurring function (i.e. the point spread function or psf). A necessary condition for the algo- 
rithm to work is that the blurring function be invariant over the entire restored image field 
(stationarity), and we note that the degradation in purely a linear operation. The general 
approach in the algoritm is to find a pair of functions whose convolution gives the input 
data within a set of physical constraints. The most important constraint requires that both 
the image and the psf be positive. When the data is noisy, this will not be strictly true, 
but minimizing the negativity in the solution appears to be sufficient in almost all cases. 
Another constraint is applied by defining a ”support” region in image space where the im- 
age (or psf) is non-zero, and then resetting that area to zero after each iteration. Finally, 
the signal-to-noise ratio in the Fourier transform (FT) of the image or psf constrains the 
dynamic range of the deconvolution operation by using a Wiener filter for the deconvolution 
operation. While it has not been proven that the restored images from BID are unique, 
complicated images appear to converge on only one sensible solution. 
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2 Implementation 


A flow diagram for the technique is given in figure 1. One starts with a degraded 
image and an initial estimate of the point spread function (psf). The initial psf can be 
randomly chosen, however the number of iterations required for the algorithm to converge is 
highly dependant on how close the first estimate of the psf is to the actual psf. Both inputs 
are Fourier transformed and a deconvolution is performed by constructing a Wiener filter 
from the FT of the psf and an estimate of the noise in the image spectrum. The technique of 
Wiener (or Optimum) filtering damps the high frequencies and minimizes the mean square 
error between each estimate and the true spectrum. 


The Weiner filtering spectrum, o/(u,v), usually takes the form: 


oj(u, v ) = o(u , v ) • 


\p(u,v)\ 2 + |n(u,v)| 2 


(1) 


o(u,v), p(u,v ) and n(u,v ) are the object, PSF and noise spectra respectively. For white 
noise, n(u, v) usually can be replaced with a constant estimated as the rms fluctuation of 
the high frequency region in the object spectrum (i.e. at frequencies where the object power 
is negligible). 

The result after Wiener filtering is transformed back to image space and positivity 
and support constraints are applied. The negatives in the image are set to zero and their 
total value summed. This value is divided by the number of pixels inside the support region 
and subtracted within that region in order to preserve the total power in the image. After 
subtraction, some areas of the image may become negative. If this is the case, the negatives 
are again truncated, summed and subtracted. This procedure is repeated until the negatives 
in the restored image are reduced to a negligible level. The original degraded image is then 
deconvolved by the restored image obtained from the first iteration. The result is transformed 
back to image space. Again, positivity and support constraints are enforced. The result is a 
new estimate of the psf. The iteration continues until a stable solution is found. A damping 
factor is used to stabilize the iteration, particularly important when the psf estimate is still 
inaccurate. About 20% of the image (or psf) from the previous cycle is averaged with the 
new image (or psf) in the early stages of the process. This percentage is reduced when the 
iteration has nearly converged. 

Two criterion have been found to be very useful in determining the completion of 
the iterations: the ratio of positive power to negative power in the restored image and psf; 
and the rms difference between images or psfs from successive iterations. Both criteria drop 
irregularly in the first few cycles of the iteration, but they both level off and stabilize when 
the operation is close to convergence. After examination of the output image and the psf, 
the results may be fed back into the loop for continued iterations. 

There are a number of parameters which must be chosen in order to ensure conver- 
gence and an optimum result. Probably the most important are accurately estimating the 
signal-to-noise ratio in the data, allowing construction of the Wiener filter, and carefully 
defining the region for the support constraint. It is also very important that the image and 
psf remain aligned with the support, since truncation seriously degrades the process. This 
is done by centering the initial image and psf, calculating the two support regions, and then 
recentering the psf after each cycle. 



Appendix A contains listings of the software developed for BID. It consists of four 
packages: the main program (bid.f), a collection of subroutines used by BID (bidsubs.f), 
a set of subroutines used to interface with the IRAF data reduction package (irafsubs.f), 
and a multi-dimensional fast fourier transform (fourt.f). Each of the packages has internal 
documentation that describes the functions of each routine and aids the selection of param- 
eters when running the program. This choice of parameters is highly dependant on the data 
that is to be deconvolved, and it must be tuned by the user through experience. However, 
suggested starting values are provided that should help initially. Also included is wiener 
deconvolution program that can be used after BID has found the psf. Reinsertion of the 
original data and the psf from BID into DCW (dcw.f) allows deconvolution of the image with 
an adjustable low-pass filter and a choice of wiener filter parameter to produce smoother or 
sharper reconstructions. Once the psf has been found, the signal-to-noise ratio in the data 
determines the frequency range and enhancement level acceptable in the deconvolved image. 
The programs must be compiled inside IRAF using the FC compliler command. 


3 Einstein X-ray Images 


We have applied BID to an extensive set of supernova remnants from the Large 
Magellenic Cloud. The results have varying success, mostly depending on the signal-to-noise 
ratios in the raw data. In the very low surface brightness SNR’s, the result is not much more 
than a quasi-optimum smoothing. All deconvolved images were produced using the same 
parameters in BID (to find the psf) and then in DCW to produce the image: we used the 
parameters suggested in the BID program and we ran BID with the psf support-only option. 
In DCW, we used a frequency cutoff of 1.0 (the maximum value) and a wiener filter multiplier 
of 4.0. Some of the images could be restored with smaller wiener multipliers (producing a 
sharper image) but we show only the results for this one value to allow comparison between 
data sets. 

Figures 1 and 2 are for Cyg X-l. The deconvolved image is point-like with the 
elongation seen in the raw data (figure 1), presumably due to aberrations removed. DEM71 
is an extended low surface brightness remnant. The deconvolved image (Figure 4) is only a 
somewhat smoothed version of the input, unsurprising for the low signal-to-noise ratio in the 
raw image (Figure 3). N23 is another extended remnant. Here the deconvolved image (Figure 
6) has sharpened and defined the bright knots, allowing a better feeling for the structural 
detail than in the raw data (Figure 5). The results for N49, N49B, N63A and N132 (Figures 
7-14) are similar to N23: some improvement in definition of the brighter features and an 
overall smoothing of the images allowing easier interpretation. 

Figures 15 and 16 show the data and deconvolution for the results obtained for N103B. 
In this case, the data are believed to suffer from image motion in the up-down direction, 
giving the image its double lobed appearance. The deconvolution has partially compensated 
for this, substantially brightening the lower ’’lobe”. However, image motion is difficult to 
completely remove without very high signal-to-noise ratios since its transfer function drops 
rapidly to zero. 

Figures 17 and 18 show the input image and deconvolution for SNR0102, a circular 
ring remnant, seen almost face-on. The deconvolution shows a number of interesting features 
including a number of bright knots and asymmetries in the ring. Particularly striking is the 



dark break in the ring near the top of the image. The raw data shows hints of this feature, 
but it is cleanly defined in the reconstruction. 

Figures 19 and 20 produce what is perhaps the most perplexing result from the 
SNR0519. This is a rather small SNR, appearing rather amorphous in the raw data, however 
the deconvolution produces a very different result: a bright central peak with a small bright 
region, all sitting on a fainter extended remnant. This data was run for several hundred 
extra cycles and with different random starting psf’s. In all cases the results wer the same. 
The image in figure 20 is after 200 iterations of BID, very nearly the same result as the 100 
iteration output. Further analysis will be required to determine whether the result is real or 
whether some special problem with data produced it. This is the only result that does not 
appear to be consistent with the original data. 

Figures 21 and 22 show the results for SNR0540, a very compact remnant. The data 
and the deconvolution show a very bright central core and some extend flux around the SNR. 
Figures 23 and 24 present two other data displays of the deconvolution, a contour plot and 
a surface plot. These displays hint at a flattening of the central core, possibly corresponding 
to the ring seen in optical images. Since it is only a few pixels wid, it would be difficult to 
conclude the exact nature of the flattening. 

These result demonstrate the capabilities and the limitation of the BID program. In 
all cases, the results appear to be stable, and in all cases but one, consistent with original 
data. For data with higher signal-to- noise ratio, substantial sharpening of the image features 
is obtained. For the very low surface brightness images, the result is close to an optimally 
smoothed result. 


4 Summary 

BID appears to be a powerful new tool for high angular resolution astronomy. While 
the technique requires fairly high signal-to-noise ratios in the data, substantial improvement 
in image sharpness can be obtained even for low flux x-ray images. The programs are 
relatively easy to using, interfacing with the IRAF package to allow reading and writing 
IRAF images. 
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Figure Captions 


Figure (1) - Einstein Data for Cyg X-l 

Figure (2) - BID Deconvolution for Cyg X-l 

Figure (3) - Einstein data for DEM71 

Figure (4) - Bid Deconvolution for DEM71 

Figure (5) - Einstein Data for N23 

Figure (6) - BID Deconvolution for N23 

Figure (7) - Einstein data for N49 

Figure (8) - Bid Deconvolution for N49 

Figure (9) - Einstein Data for N49B 

Figure (10) - BID Deconvolution for N49B 

Figure (11) - Einstein data for N63A 

Figure (12) - Bid Deconvolution for N63A 

Figure (13) - Einstein data for N103B 

Figure (14) - Bid Deconvolution for N103B 

Figure (15) - Einstein Data for N132D 

Figure (16) - BID Deconvolution for N132D 

Figure (17) - Einstein Data for SNR0102 

Figure (18) - BID Deconvolution for SNR0102 

Figure (19) - Einstein data for SNR0519 

Figure (20) - Bid Deconvolution for SNR0519 

Figure (21) - Einstein Data for SNR0540 

Figure (22) - BID Deconvolution for SNR0540 

Figure (23) - Contour Plot of the BID Deconvolution of SNR0540 

Figure (24) - Surface Plot of the BID Deconvolution of SNR0540 
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Figure (1) - Einstein Data for Cyg X-l 



ORIGINAL FAGS 

BLACK AND WHITE PHOTOGRAPH 

cygxld300w4 




Figure (2) - BID Deconvolution for Cyg X-l 
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Figure (3) - Einstein data for DEM71 
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Figure (4) - Bid Deconvolution for DEM71 
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Figure (5) - Einstein Data for N23 
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Figure (6) - BID Deconvolution for N23 
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Figure (7) - Einstein data for N49 
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Figure (8) - Bid Deconvolution for N49 
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Figure (9) - Einstein Data for N49B 
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Figure (10) - BID Deconvolution for N49B 
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Figure (11) - Einstein data for N63A 
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Figure (12) - Bid Deconvolution for N63A 
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Figure (13) - Einstein data for N103B 
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Figure (14) - Bid Deconvolution for N103B 
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Figure (15) - Einstein Data for N132D 
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Figure (16) - BID Deconvolution for N132D 




snr0102gzm 


ORIGINAL PAGE 

BLACK AND WHITE PHOTOGRAPH 



Figure (17) - Einstein Data for SNR0102 
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Figure (18) - BID Deconvolution for SNR0102 
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Figure (19) - Einstein data for SNR0519 
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Figure (20) - Bid Deconvolution for SNR0519 
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Figure (21) - Einstein Data for SNR0540 
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Figure (22) - BID Deconvolution for SNR0540 
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Figure (23) - Contour Plot of the BID Deconvolution of SNR0540 
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Figure (24) - Surface Plot of the BID Deconvolution of SNR0540 
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Appendix A 



program bid 

c last edited: January 6, 1992 
c Compile in IRAF using FC 
c link with bidsubs, irafsubs, and fourt 
c perform blind iterative deconvolution 
c using wiener filters only 
c needs input image, guess at psf 
c 

c BID is set up to take square IRAF image arrays of up to 512x512 
c Smaller arrays or smaller processing areas are chosen by 
c entering the parameter KSZ. 

c Images smaller than KSZ are set in KSZ arrays and the outer 
c regions are filled with either zero or and average constant, 
c Images larger than KSZ (but smaller than 513x513) and 
c chunked to KSZ area, at an entered center, 
c File names are for IRAF files without the .imh. 
c 1st gue3S psf can be either a file or it can be generated with 
c PSFGEN . 

c Program runs for MAXC cycles. 

c BETA is the percentage (0 to 1.0) of the previous cycle averaged 
c with the current cycle. Usually 0.8 near the start of the 
c process, but changed to 1.0 near convergence, 
c Resulting psf's and images can be read back into the program 
c to perform additional cycles. 

c CF is the radius in frequency space used for the Wiener filter 
c noise estimate. 1.0 is a good starting point, 
c W1 is the Wiener filter scale factor. 1.0 is a good starting 
c value. Program DCW allows generation of Wiener filtered 
c versions of the image using the psf from BID, but with different 
c values for wi . Larger values for W1 give a smoother image, 
c smaller values enhance the high frequencies, 
c Output file names should be given without the trailing .imh. 
c Object and psf masks are generated using program MAKEMASK. 
c Recommended starting masks should be of order twice the diameter 
c of the data and expected psf diameters. For an object that nearly 
c fills the field, a constant mask (equivalent to no mask) should 
c be used for the object and only a psf mask is used, 
c 

c If you encounter problems in using this program, contact 
c Peter Nisenson, Center for Astrophysics, (617) 495-7394 
c 


parameter (isz=512, iszl=isz/2, isz2=iszl+l) 

c complex fourier space arrays : c = g*f 
complex c(isz*isz) 
complex g( isz* isz) 
complex f(isz*isz) 

c real image space arrays, 1-d here but all routines think they are 2- 
real cob j (isz*isz) , gref (isz*isz) , fout (isz*isz) , last f out (isz 
c masks and a temporary buffer 

real maskr(isz*isz) ,mask(isz*isz) , tempb (isz*isz) 

c total power of input image 
real power 

c two working arrays, one for fourt, one for loadim 
real work (2*isz) , rwk (isz) 
c testo, testr, rmsdif f array, for output to file 
real stats (isz, 3) 

c temp integer array for loadim : 

integer iph(isz,isz) 


T) * 



c nzero and nzeror are number of non-zeros (ones) in mask and maskr : 
integer nzero, nzeror 

c parameters for fourt : 

integer*4 nn (2) , ndim, isign, if orm 

c file names : 

character namel*80, name2*80, name3*80, name4*80, name5*80 

print *, ' BID : iterative deconvolution given initial ' 
print *, ' guess at the psf, uses wiener deconvolution, ' 
print *, ' support and positivity constraints to iteratively' 
print *, ' approximate the correct psf . ..' 
print *, ' ' 

print *, ' Enter log file name' 
read (*,*) name5 

open (unit=l, name=name5, type=' new' , f orm=' formatted' ) 
c initialize . . . 

print *, ' enter array dimension (ksz) ' 
read (*,*) ksz 

write (1,*) ' array dimension ',ksz 

kszl=ksz/2 
ksz2=ksz*ksz 
nn (1) =ksz 
nn (2) =ksz 
ndim=2 
c 

c read in image file, save total power for rescaling output 

c 

print *, ' enter input image file name ' 
read (*,'(a)') namel 
write (1,*) ' Input file name', namel 
call loadim (namel, cob j , iph, rwk, ksz, 3) 

print *, ' 1 for border square mask, 2 for circular, 0 for none' 
read (*,*) isqm 
if ( isqm. eq. 1) then 

print *, 'Enter the sqmask border width, 0 for none' 
read (*,*) iradm 
elseif (isqm.eq.2) then 

print *, ' Enter circular mask radius ' 
read (*,*) radm 
endif 

c apply square gaussian border 

if (isqm. eq. 1) call sqmask (cob j, iradm, ksz) 

c apply circular gaussian border 

if (isqm.eq.2) call circmask (cob j, radm, ksz) 


power=0 . 0 
do j=l,ksz2 

power=power+abs (cobj ( j) ) 
enddo 

print *, ' Total power in image :', power 
write (1,*) ' Total power in image ', power 
c 

c read in the psf file and center it, scale total power to 1.0 
c 

print *, ' 1 to enter psf file, 2 to generate gaussian' 
read(*,*) igss 
if (igss .eq. 1) then 

print *, ' enter input psf file name' 



read (*,'(a)') namel 

write (1,*) 'Input psf file name namel 
call loadim (namel, gref, iph, rwk, ksz, 3) 
call centre (gref , tempb, ksz, kszl) 
call scale (gref , 1 . 0, ksz) 
else 

c generate 1st pass psf 
c 

print *, 'Generate psf' 

print *, ' Enter 1/e radius for psf ' 

read (*,*) wid 

write (1,*) ' Gaussian psf 1/e radius ',wid 
krn=l 

cf=l . 0*kszl 

call gauss (gref , wid, krn, ksz) 

if orm=0 

isign=-l 

call copyr (g, gref , isign, ksz) 

call f ourt (g, nn, ndim, isign, if orm, work) 

call rollof f (g, cf , ksz) 

if orm=l 

isign=l 

call f ourt (g, nn, ndim, isign, if orm, work) 
call copyr (g, gref , isign, ksz) 

call f fmax (amax, amin, imax, jmax, imin, jmin, gref, ksz) 

amm=amax-amin 

do j=l,ksz2 

gref ( j) = (gref ( j) -amin) /amm 
enddo 
endif 
c 

c enter reconstruction parameters 
c 

print *, ' enter max # cycles for iteration ' 
read (*,*) maxc 

write (1,*) ' max # iterations ' ,maxc 

print *, ' enter beta, the damping factor (try 0.8)' 

read (*,*) beta 

write (1,*) ' beta ',beta 

print *, ' enter radius for noise estimate (0 to 1.0)' 
read (*,*) cf 

write (1,*) ' noise radius ',cf 
cf = cf*kszl 
icf = cf 


print *, ' enter object wiener filter parameter (Try 1.0)' 
read (*,*) wol 

write (1,*) ' object Wiener filter multiplier ',wol 

print *, ' enter psf wiener filter parameter (Try 1.0)' 
read (*,*) wpl 

write (1,*) ' psf Wiener filter multiplier ',wpl 

print *, '1 for circular object support, 2 for threshold, 

+ 3 for none' 

read (*,*) ispp 

suppo=0 . 
thresh=0 . 0 



if (ispp .eq. 1) then 

print *,' enter the object support radius, 0 for none' 
read (*,*) suppo 

write (1,*) ' object support radius suppo 
endif 

if (ispp .eq. 2) then 

print *, ' enter threshold level for support ' 
read (*,*) thresh 

write (1,*) ' object support threshold thresh 
endif 

print ' 1 for psf positivity, 2 for just support' 
read (*,*) ipss 

print *,' Enter the psf circular support radius, 0 for none' 
read (*,*) suppr 

write (1,*) ' psf support radius ', suppr 

print *, ' enter output image file name ' 
read (*,'(a)') name2 

write (1,*) ' enter output image file name ',name2 

print *, ' enter output psf file name ' 
read (*,'(a)') name3 

write (1,*) ' output psf file name ',name3 

print *, ' enter output statistics file name ' 
read (*,'(a)') name4 

write (1,*) ' Output stat file name ',name4 
c 

c load previous output image or set f,ff to 0.0 
c 

print *, ' 1 to enter output image from previous iteration' 
read (*,*) iobj 
if (iobj .eq. 1) then 

print *, ' enter file name' 
read (*,'(a)') namel 

write (1,*) ' Old image file name ', namel 
call loadim(namel, fout, iph, rwk, ksz, 3) 
isign=-l 
if orm=0 

call copyr (f, fout, isign, ksz) 

call fourt (f , nn, ndim, isign, if orm, work) 

call normf(f,ksz) 

c ft back to image space so lastfout will be scaled correctly 
c for rmsdiff calculations . . . 
isign=l 
iform=l 

call fourt (f , nn, ndim, isign, iform, work) 
call copyr (f, lastfout, isign, ksz) 

isign=-l 

iform=0 

call fourt (f , nn, ndim, isign, iform, work) 
call normf(f,ksz) 
else 

call zero(f,ksz) 
do j=l,ksz2 

lastfout ( j) =0 . 0 
enddo 
endif 


c 


print *, ' generate masks' 
if (ispp .eq. 1 .and. suppo .ne . 0 . ) then 



call maskg (mask, suppo, ksz) 
elseif ( ispp . eq . 2 ) then 
do j=l,ksz2 

if (cobj(j) .gt .thresh) then 
mask(j) = 1.0 

else 

mask ( j ) =0 . 0 
endif 
enddo 
else 

do j=l,ksz2 
mask(j) = 1.0 
enddo 

endif 

c setting psf mask to radius of "supp" 
if (suppr .ne . 0 . ) then 

call maskg (maskr, suppr, ksz) 

else 

do j=l,ksz2 

maskr (j) = 1.0 
enddo 

endif 


c now count number of l's in masks 
nzero = 0 
nzeror= 0 
do j=l,ksz2 

if (mask ( j ) . ne . 0 . ) nzero=nzero+l 
if (maskr ( j ) . ne . 0 . ) nzeror=nzeror+l 
enddo 

print *, ' nzero, nzeror ' , nzero, nzeror 
c 

c ft and normalize image for first deconvolution pass 
c 

if orm=0 
isign=-l 

call copyr (c, cob j , isign, ksz) 
call f ourt (c, nn, ndim, isign, iform, work) 
call normf (c,ksz) 
c 

c perform deconvolution in loop 
c 

print *, ' deconvolving' 

do 1000 iz=l,maxc 
c ft, shift and normalize psf . . . 
iform=0 
isign=-l 

call copyr (g, gref, isign, ksz) 

call f ourt (g, nn, ndim, isign, iform, work) 

call shift (g, ksz) 

call normf (g, ksz) 

c output = input deconvolved by psf : 
c if first time through, deconvolve with beta=1.0 
if (iob j . ne . 1 . and. iz .eq. 1) then 

call wienerf (c, g, f , ksz, cf , wol, 1 . 0) 
else 

call wienerf (c,g, f , ksz, cf, wol, beta) 
endif 

c ft output to image space, apply positivity, ft back, normalize 
if orm=l 
isign=l 

call f ourt (f , nn, ndim, isign, iform, work) 



call copyr (f, fout, isign, ksz) 

call posf (fout, mask, nzero,testo,testno, ksz) 
print *, ' ' 

print *, ' Iteration #',iz 

write (1,*) ' ' 

write (1,*) 'Iteration # ',iz 

print *, 'obj: summed neg/pos ',testo,' #neg/#pos ',testno 
write (1,*) 'obj: sum neg/pos ',testo, ' #neg/#pos ' , testno 

c psf = input deconvolved by output 

iform=0 

isign=-l 

call copyr (f, fout, isign, ksz) 

call fourt (f , nn, ndim, isign, iform, work) 

call normf (f , ksz) 

call wiener f (c, f , g, ksz, cf, wpl, beta) 

call shift (g, ksz) 

iform=l 

isign=l 

c ft psf to image space, centre, apply positivity 
call fourt (g, nn, ndim, isign, iform, work) 
call copyr (g, gref, isign, ksz) 
call centre (gref , tempb, ksz, kszl) 

if (ipss .eq. 1) then 

call posf (gref ,maskr, nzeror, testr, testnr, ksz) 
else 

call support (gref,maskr, testr, testnr, ksz2) 
endif 

print *, 'psf: summed neg/pos ', testr,' #negs/#pos ', testnr 
write (1,*) 'psf: sum neg/pos ', testr,' #negs/#pos ', testnr 
c 

c testo is sum of negatives over sum of positives in output 
c prior to positivity, testr is same for psf . . . 
c 

rmd = rmsdiff (lastfout, fout, ksz) 
print *, ' rmsdiff : ' , rmd 
write (1,*) ' rmsdiff : ' , rmd 
stats (iz, 1) =testo 
stats (iz, 2) =testr 
stats (iz, 3) =rmd 

1000 continue 

c 

c stop : scale total power of output to input power, psf to 1.0 
c then save 
c 

2000 continue 

call scale (fout, power, ksz) 
rtot=0 . 0 
do j=l,ksz2 

rtot=rtot+fout ( j) 
enddo 

print *, ' total power in image ' , rtot 
call saveim(name2, fout, ksz, 3) 

call scale (gref , 1 . 0, ksz) 
call saveim(name3,gref , ksz, 3) 

call saverect (name4, stats, ksz, 3, ksz, 3) 


stop 

end 



program dew 

c last edited: November 4, 1991 
c Compile in IRAF using FC 
c link with: bidsubs, irafsubs, and fourt 

c uses wiener deconvolution 

c Allows application of a edge mask for images which extend 
c beyond the frame boundary 

c Best results usually obtained with square low-pass filter 
c and cutoff of 1.0 

c Smoothness of output image is adjusted by chnaging the wiener 
c multiplier - greater value for smoother image, 
c 
c 

parameter (isz=512, iszl=isz/2, isz2=iszl+l) 

real*4 work(2*isz) 

real*4 ph (isz*isz) , rwk (isz) 

complex sne (isz*isz) , refc (isz*isz) , sneft (isz*isz) , snede (isz*isz) 
integer nn(2) 
integer iph(isz*isz) 

character namel*80, name2*80, name3*80 


print * 
print * 
print * 
print * 
print * 
print * 
print * 


deconvolution using wiener filter ' 

wiener parameter wl : 0 => straight deconvolution 

wl >> 0 gives noise reduction ' 

can truncate negatives at end or apply the ' 

fienup algorithm to the output ' 


print * 
read (* 
print * 
read < * 


ksz2=ksz*ksz 


' enter input object file name 
' (a ) ' ) namel 

' enter image dimension ' 

*) ksz 


9 


call loadim (namel, ph, iph, rwk, ksz, 3) 

print *, ' 1 for border square mask, 2 for circular, 0 for none' 
read (*,*) isqm 
if ( isqm. eq. 1) then 

print *, 'Enter the sqmask border width, 0 for none' 
read (*,*) iradm 
elseif (isqm.eq.2) then 

print *, ' Enter circular mask radius ' 
read (*,*) radm 
endif 

c apply square gaussian border 

if (isqm. eq. 1) call sqmask (ph, iradm, ksz) 

c apply circular gaussian border 

if (isqm. eq. 2) call circmask (ph, radm, ksz) 


rintot=0 . 0 
do j=l,ksz2 

sne ( j ) =ph ( j ) 
rintot=rintot+ph ( j ) 

enddo 

print *, ' total power in image ' , rintot 

print *, ' enter input reference file name ' 

read (*, ' (a) ' ) name2 

call loadim (name2 , ph, iph, rwk, ksz, 3) 

do j=l,ksz2 



enddo 


ref c ( j) =ph ( j) 


print *, ' enter wiener filter multiplier, 1 is default' 
read (*, *) wl 

print *, ' enter radius for wiener noise estimate (0 to 1.0) ' 
read (*,*) wr 
wr=wr*ksz/2 . 0 

print *, ' enter low frequency cutoff, (0 to 1.0) ' 
read (*,*) cf 
cf=cf *ksz/2 . 0 

print *, ' 1 for circular low pass filter, 2 for square filter, 
read (*,*) iflt 

c set up parameters for ft 
ndim=2 
iform=0 
nn (1) =ksz 
nn (2) =ksz 
isign=-l 

c ft image 

print *, ' transforming image . ..' 
call f ourt (snc, nn, ndim, isign, if orm, work) 

c ft psf 

print *, ' transforming reference . ..' 
call fourt (refc, nn, ndim, isign, iform, work) 

do j=l,ksz2 

sncft ( j)=snc ( j) 

enddo 

c perform wiener deconvolution 

600 print *, ' doing wiener deconvolution . ..' 

call wienerd (snc, refc, ksz, wr, wl) 
do j=l,ksz2 

sncdc( j)=snc(j) 

enddo 

c low pass filter 

650 if (cf ,ne . 0 . ) call f ilt (snc, ksz, cf , if It) 

c shift image to center of field 
call shift (snc, ksz) 

if orm=l 
isign=l 

c ft back to image space 

print *, ' transforming back to image space ...' 

call fourt (snc, nn, ndim, isign, iform, work) 

do j=l,ksz2 

if (real (snc (j)) . lt.0.0)snc(j)=0.0 
enddo 

do j=l,ksz2 

ph ( j ) =real (snc ( j) ) 
enddo 

routtot=0 . 0 
do j=l,ksz2 

routtot=routtot+ph ( j ) 

enddo 

print *, 'output power before renorm ' , routtot 


or 0' 



rsf act=rintot/ rout tot 
do j=l,ksz2 

ph { j ) =ph ( j) * rsf act 

enddo 

routtot=0 . 0 
do j=l,ksz2 

routtot=routtot+ph ( j ) 

enddo 

print *, ' total output power ' , routout 
print *, ' enter output file name ' 
read (*,'(a)') name3 
call saveim(name3,ph, ksz, 3) 


c again? 

print *, ' 0 to quit, 1 to redo with new wl, 2 to redo with new cf 
read (*,*) ilp 
if (ilp .eq. 0) goto 99999 
if (ilp .eq. 1) then 

print *, ' enter wiener filter multiplier, 1 is default' 
read (*,*) wl 

print *, ' enter radius for wiener noise estimate (0 to 1.0) ' 
read (*,*) wr 
wr = wr*ksz/2.0 

print *, ' enter low frequency cutoff, (0 to 1.0)' 
read (*,*) cf 
cf = cf*ksz/2.0 

print *, ' 1 for circular low pass filter, 2 for square filter 
read (*,*) if It 
do j=l,ksz2 

snc ( j) =sncft ( j) 

enddo 
goto 600 
endif 

if (ilp .eq. 2) then 

print *, ' enter low frequency cutoff (0 to 1.0)' 
read (*,*) cf 
cf = cf*ksz/2.0 

print *, ' 1 for circular low pass filter, 2 for square filter 
read (*,*) iflt 
do j=l,ksz2 

snc ( j) =sncdc ( j) 

enddo 
goto 650 
endif 


99999 stop 
end 



c Irafsubs.f 

c loadim : load an image (iraf, real or integer) 

c saveim : save an image (iraf or real) 

c saverect : save a rectangle chunked from corner of larger 

c rectangle to iraf file (for saving convergence statistics) 

subroutine loadim (name, a, itmp, rwk, isz, ftyp) 

c last edited july30,1990 
c loads iraf files into fortran program 

c itmp is an integer array, iszxisz, for loading integer 
c images, rwk is a one-row working real array, 

c ftyp is file-type : 1 means integer, 2 means real, 3 means iraf. 
c 

character name* 80 
real a (isz, isz) 
integer itmp (isz, isz) 
real rwk (isz) 
integer ftyp 
integer im 
integer axl(7) 
integer naxis 
integer dtype 
integer ier 

integer xcentr, ycentr, il, i2, jl, j2 
integer ist, iend, jst, jend, idim, jdim 
real bsum 
isz4=isz*isz*2 
isz6=isz4*2 
c 
c 

if (ftyp .eq. 1) then 

open (unit-1, f ile=name, access=' direct' , status=' old' , 

+ form=' unformatted' , recl=isz6) 

read (unit=l, rec=l) itmp 

open(unit=l, file=name, status=' old' , form=' unformatted' , 

+ access=' direct' , recl=isz6) 

read(unit=l) itmp 
do j=l,isz 
do i=l,isz 

a (i, j ) = itmp ( i , j) 
end do 
end do 

close (unit=l) ( 

return 

else if (ftyp .eq. 2) then 

open (unit=l, f ile=name, access=' direct' , status=' old' , 

+ form=' unformatted' , recl=isz6) 

read (unit=l, rec=l) a 
close (unit=l) 
return 

else if (ftyp .eq. 3) then 

call imopen (name, 1, im, ier) 
if (ier .ne. 0) goto 99999 
call imgsiz (im, axl, naxis, dtype, ier) 
if (ier .ne. 0) goto 99999 

if ((axl(l) .gt. isz) .or. (axl(2) .gt. isz)) then 
print *, ' image bigger than array size :' 
print *, ' axl(l): ', axl(l),' ... axl (2): ', axl (2) 
write (*,*) 'input centre of extraction box' 
read (*,*) xcentr, ycentr 

xcentr = min ((axl(l) - isz/2) , max (isz/2, xcentr)) 



& 


500 


ycentr = min ((axl (2) - isz/2) , max (isz/2, ycentr) ) 

11 = xcentr - isz / 2 + 1 

12 = il + isz - 1 

jl = ycentr - isz / 2 + 1 
j2 = jl + isz - 1 

call imgs2r (im, a, il, i2, jl, j2, ier) 
if (ier .ne. 0) goto 99999 

else if ((axl(l) .It. isz) .or. (axl(2) .It. isz)) then 
idim = axl (1) 
jdim = axl (2) 
ist = (isz-idim) /2 + 1 
iend = ist+axl(l)-l 
jst = (isz- jdim) /2 + 1 
jend = jst+axl(2)-l 
do j=jst, jend 

call imgs2r (im, rwk, 1, axl (1) , 1+ j- jst, 1+ j- jst, ier) 
if (ier .ne. 0) goto 99999 
do i=ist,iend 

a (i, j) =rwk (i-ist+1) 
end do 
end do 
bsum=0 . 0 

print *, ' image smaller than array size 
print *, 

1 to calculate average background, 0 to set to zero' 
read (*,*) iback 
if (iback .eq. 0) goto 500 
do i=ist,iend 

bsum = bsum + a (i, jst) +a (i, jend) 
end do 

do j=jst,jend 

bsum = bsum + a (ist, j) +a (iend, j) 
end do 

bsum = bsum/ (2 .* jdim + 2.*idim) 

write (*,*) ' average background : ',bsum 

if (jst .gt. 1) then 
do j=l,jst-l 
do i=l,isz 

a (i, j ) =bsum 
end do 
end do 
end if 

if (ist .gt. 1) then 
do i=l,ist-l 

do j=jst, jend 
a (i, j ) =bsum 
end do 
end do 
end if 

if (iend .It. isz) then 
do i=iend+l,isz 
do j=jst,jend 

a (i, j) =bsum 
end do 
end do 
end if 

if (jend .It. isz) then 
do j=jend+l,isz 
do i=l,isz 

a (i, j) =bsum 
end do 
end do 
end if 



else 

11 = 1 

12 = axl(l) 
jl = 1 

j2 = axl (2) 

call imgs2r (im, a, il, i2, jl, j2, ier) 

to 99999 

end if 

call imclos (im, ier) 
if (ier .ne. 0) goto 99999 
return 
end if 
c 

99999 call imemsg (ier, name) 
write (*, 2222) name 
2222 format ('error: ',a80) 

stop 
end 
c 

c 

c 

subroutine saveim(name, a, isz, ftyp) 


if (ier .ne. 0) go 


c last edited: July 30, 1991 
c saves images as iraf files 

character name*80 
real a(isz,isz) 
integer isz 
integer ftyp 

integer axl (7) ,naxis, im, ier 
if (ftyp .eq. 3) then 
axl (1) =isz 
axl (2) =isz 
naxis=2 

call imcrea (name, axl, naxis, 6, ier) 
if (ier .ne. 0) goto 99998 
call imopen (name, 3, im, ier) 
if (ier .ne. 0) goto 99998 
call imps2r (im, a, 1, isz, 1, isz, ier) 
if (ier .ne. 0) goto 99998 
call imclos (im, ier) 
if (ier .ne. 0) goto 99998 
return 

open (unit=l, f ile=name, status=' new' , form=' unformatted' , 
access=' direct' , recl=4*isz*isz) 
write (unit=l, rec=l) a 
close (unit=l) 
return 

99998 call imemsg (ier, name) 
write (*, 2222) name 
2222 format ('error: ',a80) 

stop 
end 


else 

+ 


end if 


Q ******************************************************************* 

c 

subroutine saverect (name,a,dl,d2, rdl, rd2) 

c save a few lines of data for convergence plots 
character name*80 
integer dl, d2, rdl, rd2 



99997 


c 


real a(rdl,rd2) 

integer axl (7) , naxis, im, ier 

print *, ' dl : ' , dl, ' , d2:',d2,', rdl:', rdl,', rd2 : ' , rd2 
axl (1) =dl 
axl (2)=d2 
naxis=2 

call imcrea (name, axl, naxis, 6, ier) 

if (ier .ne. 0) goto 99997 

call imopen (name, 3, im, ier) 

if (ier .ne. 0) goto 99997 

call imps2r (im, a, l,dl, l,d2, ier) 

if (ier .ne. 0) goto 99997 

call imclos <im,ier) 

if (ier .ne. 0) goto 99997 

return 

call imemsg (ier, name) 

write (*,' ("error: ", a80)') name 

stop 

end 



c 

c bidsubs.f : subroutines for idconf and makemask, etc... 
c 

c last edited December 6, 1991 
c 

c 

c routines are : 

c cent : Finds center of mass 

c centre : centres the object using ffmax to find psf center 

c circmask: multiply image by circular rolloff mask 

c copyr : copies real array into complex array or vice versa, 

c depending on value of isign 

c cwind : a rolloff window function in f-space 

c dtr: removes trend from solar images before BID 

c ffmax : finds position and values of minimum and maximum pixels 

c in image 

c filt : apply square or circular rolloff filter using cwind at 

c given radius . . . 

c maskg: Generates gaussian masked image 

c normf : normalizes in fourier space so a(l,l)=l 

c pcent : print array center 

c posf : the old version (and best working, it seems) of pos, 

c iteratively truncates negatives and adds them back in ... 

c rmsdiff : rmsdiff between two images (real) , then stuffs second 

c image into first . . . 

c rolloff : rolls off high frequencies using cwind 

c scale : scales image to given power level 

c shift : reverses sign of every other element in ft, shifts 

c object from center to corners or vice versa . . . 

c sqmask: Generate edge rolloff mask for solar images 

c subtract: Subtracts low frequencies for solar detrending 

c support: Applies image plane support constraint for BID 

c wienerd : wiener filter deconvolution, with averaging in of 

c previous f (f = ( 1 . 0-beta) *f+beta*c/g) 

c wienerf : another version of wienerd 

c zero : zeros a complex array 

c 

Q ***************★*★★★**★★*★**★*★****★*★*★★★★**★*★★**★★*★******★*★***** 

C 

subroutine cent (ci, c j , tm, ph, isz) 
c 

c calculates center of mass of iszxisz image 
c 

c xm = x coordinate of centre of mass 
c ym = y coordinate of centre of mass 
c tm = running sum of values in image array 
c ph() = image array 

c isz = characteristic image dimension (ie. 128) 
c 

dimension ph (isz, isz) 
real copy 

call ffmax (amax, amin, imax, jmax, imin, jmin, ph, 
amm=0 . 3*amax 

zl=isz/2 
z2=iszl+l 
tm=0 . 
xcl=0 . 
ycl=0 . 
do j=l,isz 
do i=l,isz 

copy = ph (i, j) 


isz) 



if ( copy . ge . amm) then 
tm=tm+copy 
xcl=copy*i+xcl 
ycl=copy* j+ycl 
endif 
enddo 
enddo 

ci=xcl/tm 
c j=ycl/tm 

print *, ' center at ',ci,cj 

return 

end 

c 

Q ★★*★★★★★★★★★★★***★★★★★*★*★**★**★*★★**★★★★*★★**★★*★★'**★****★★★***★★ 

C 

subroutine centre (a, b, isz, iszl) 
real a(isz, isz), b(isz,isz) 

call cent (ci, cj, t, a, isz) 

do j = 1, isz 

do i = 1, isz 

b (i, j ) = 0.0 
end do 
end do 

jr = (nint(cj) - iszl) - .5 
ir = (nint(ci) - iszl) - .5 
do j = 1, isz 

do i = 1, isz 

il = i + ir 

c centre object 

jl = j + jr 

if ( ( ( (il .le. isz) .and. (il .ge. 1)) .and. (jl . le . isz)) .and. 
& ( jl .ge. 1) ) then 
b (i, j) = a (il, jl) 
end if 
end do 
end do 

do j = 1, isz 

do i = 1, isz 

a (i, j) = b(i, j) 
end do 
end do 
return 
end 
c 

subroutine circmask (ob j , radm, ksz) 
real obj (ksz, ksz) , radm 

kszl=ksz/2 
ksz2=kszl+l 
rr=ksz2-radm 
if (rr . It . 1) rr=l . 
do j=l,ksz 
do i=l,ksz 

r=sqrt (1.* (ksz2-i) **2+1 . * (ksz2- j) **2) 
if ( r . ge . radm) then 

obj (i, j) =obj (i, j) *cwind(abs (r-radm) /rr) 



endif 


enddo 

enddo 

return 

end 

Q ******************************************************** 


subroutine copyr (a, ar, isign, isz) 
complex a (isz, isz) 
real ar (isz, isz) 
if (isign .eq. 1) then 
c copy into real array for pos 
do j=l,isz 
do i=l,isz 

ar (i, j) =real (a (i, j) ) 
enddo 
enddo 
else 

c copy into real array for pos 
do j=l,isz 
do i=l,isz 

a (i, j) =ar (i, j) 
enddo 
enddo 
endif 

return 

end 

c 

c ********************************************************************* 
function cwind(relr) 

c A window for fft's that has low sidelobes 

real c(4) 

data c / .074, .302, .233, .390 / 
if (relr .It. 0.) stop ' cwind' 
c careful . . . 

if (relr .gt. 1.5) then 

cwind = 0 . 

else 

c (l-r**2) , r=0 . . . 1.0 

r = 1. - (relr * relr) 
r2 = r * r 
r3 = r2 * r 

cwind = ( (c (1) + (r * c(2))) + (r2 * c{3))) + (r3 * c(4)) 
c let go all the way to zero 

if (cwind .It. 0.) cwind = 0. 
end if 
return 
end 
c 

c***** ***************** *************************** ************* 


subroutine dtr (a, work, ac, acf , ksz, fc, dc) 
c subtract out the low frequency trend in an image, then add dc to 
c make image positive . . . 

real a (ksz*ksz) , work (2*ksz) 
complex ac (ksz*ksz) , acf (ksz*ksz) 
integer*4 nn (2) , ndim, isign, iform 



power=0 . 0 
ksz2=ksz*ksz 
do j=l,ksz2 

ac ( jj = (0 . , 0 . ) 
ac ( j) =a ( j) 
enddo 
c 

ndim=2 
nn (1) =ksz 
nn (2) =ksz 
iform=0 
isign=-l 

call f ourt (ac, nn, ndim, isign, iform, work) 
do j=l,ksz2 

acf ( j)=ac ( j) 
enddo 

call f ilt (acf , ksz, fc, 2) 
call subtract (ac, acf , ksz) 
iform=l 
isign=l 

call fourt (ac, nn, ndim, isign, iform, work) 
do j=l,ksz2 

a ( j)=real (ac ( j) /ksz2) 
enddo 

call ffmax (amax, amin, imax, jmax, imin, jmin, a, ksz) 
do j=l,ksz2 

a ( j) =a ( j) -amin 

enddo 

dc=-amin 

return 

end 

c 

q *************************************************************** 

c 

subroutine ffmax (amax, amin, imax, jmax, imin, jmin, ph, isz) 

dimension ph(isz, isz) 

amax = ph (1,1) 

amin = ph (1, 1) 

do j = 1, isz 

do i = 1, isz 

if (ph(i,j) .gt. amax) then 

amax = ph (i, j ) 

imax = i 

jmax = j 

end if 

if (ph(i,j) .It. amin) then 

amin = ph (i, j ) 

imin = i 

jmin = j 

end if 

end do 

end do 

return 

end 


c 

Q ★★★★★★★★★★'it*********************************'****'****** 

c 

subroutine f ilt (cob j, isz, cf , isq) 
complex cob j (isz, isz) 
iszl=isz/2 

c low pass filter 
c 



650 if (cf .ne . 0 . 0 . and. isq.eq. 1) then 

c circular filter 

print *, 'circular low pass filter 
do 100 j=l,iszl 
do 100 i=l,iszl 

dl=sqrt ( (i-1 . ) **2+ ( j-1 . ) **2) /cf 
d2=sqrt { (i-1 . ) **2+ j**2) /cf 
d3=sqrt (i**2+( j-1. ) **2) /cf 
d4=sqrt (l.*i**2+j**2) /cf 
cobj (i, j) =cwind(dl) *cobj (i, j) 
cobj (i, isz+l-j) =cwind (d2) *cob j (i, isz+1- j ) 
cobj (isz+l-i, j) =cwind(d3) *cob j (isz+l-i, j) 
cobj (isz+l-i, isz+1- j) =cwind(d4) *cobj (isz+l-i, isz+l-j) 
100 continue 

endif 
c 

if (cf .ne . 0 . 0 . and. isq. eq. 2) then 
c square filter 

print *, 'square low pass filter' 
do 200 j=l,iszl 
do 200 i=l,iszl 

dlx= (i-1 . ) /cf 
dly= (j-1 . ) /cf 
d2x= (i-1 . ) /cf 
d2y= j/cf 
d3x=i/cf 
d3y= ( j-1 . ) /cf 
d4x=i/cf 
d4y= j/cf 

cobj (i, j) =cwind(dlx) *cwind(dly) *cobj (i, j) 
cobj (i, isz+l-j) =cwind(d2x) *cwind(d2y) *cobj (i, isz+l-j) 
cobj (isz+l-i, j ) =cwind (d3x) *cwind (d3y) *cob j (isz+l-i, j ) 
cobj (isz+l-i, isz+l-j) =cwind (d4x) *cwind (d4y) 

+ *cobj (isz+l-i, isz+1- j) 

200 continue 

endif 
c 
c 


return 

end 


c 


subroutine gauss (a, wid, krn, isz) 
real a (isz, isz) 
real wid 
integer seed 
seed=13972 
iszl=isz/2 
isz2=iszl+l 


if (krn .eq. 1) then 
x=rand(seed) 
do j=l, 97 
x=rand(0) 
enddo 
endif 


wid2=wid*wid 


do j=l,isz 
do i=l,isz 



r2= ( (isz2-i) **2+{isz2-j) **2) 
r3=r2/wid2 

if ( r3 .le. 50) then 
if(krn.eq.l) then 
rr=rand(0) 
else 
rr=l . 
endif 

a (i, j) =exp (-r2/wid) *rr 
endif 
enddo 
enddo 
return 
end 

subroutine maskg (m, rd, ksz) 
real m(ksz,ksz) 

kszl=ksz/2 
do j=l,ksz 
do i=l,ksz 

r=sqrt (float (kszl-i) **2+float (kszl- j) **2) 
if ( r . le . rd) then 
m(i, j) =1 . 0 
else 

m(i, j ) =0 . 0 
endif 

enddo 

enddo 

return 

end 

Q ★**********★★**★★★**★★*★★★*★**★***★*★★*★**★★***★★****★** 

c 

subroutine normf (a, isz) 
complex a (isz, isz) 
al = cabs (a (1, 1) ) 
if (al .eq. 0.0) then 

print *, ' error in normf: a (1,1) =0.0' 
return 
endif 

do j = 1, isz 
do i = 1, isz 

a (i, j) = a (i, j) / al 
end do 
end do 
return 
end 

c 

Q******************************************************** 

subroutine pcent(a,ksz) 
real a(ksz,ksz) 
do j=ksz/2-3, ksz/2+4 

print 111, (a (i, j) , i=ksz/2-3, ksz/2+4) 
enddo 

111 f ormat ( 8el0 . 2) 

return 
end 


Q fc******************************************************************** 

c 

subroutine posf(a, b, nzero, test, testn, isz) 



C modify 12/23/1991 


real a(isz, isz) , b(isz, isz) 
integer nzero 

c 1st apply support constraint 
if (nzero .ne . 0) then 
do j = 1, isz 
do i = 1, isz 

a (i, j) = (a (i, j) * b(i,j)) / nzero 
end do 
end do 
else 

print *, ' mask is all zeros' 
return 
endif 

c then truncate negatives 
do jz=l,3 

sump = 0.0 
sumn = 0.0 
npos = 0 
nneg = 0 
do j = 1, isz 
do i = 1, isz 

if (a(i,j) .It. 0.0) then 
sumn = sumn + a(i,j) 
a (i, j) = 0.0 
nneg=nneg+l 

elseif (a(i,j) .gt. 0.0) then 
sump = sump + a ( i , j ) 
npos = npos + 1 
end if 
end do 
end do 

if (jz .eq. l)then 
if (sump.ne . 0) then 
test = abs(sumn / sump) 
else 

print Image all negs' 

return 
endif 

if (npos . ne . 0 ) then 
testn=f loat (nneg) /npos 
else 

print *, ' No positives in image' 
endif 
endif 

c now add in negatives to keep energy constant 
ss = sumn /nzero 
do j = 1, isz 
do i = 1, isz 

if (b(i,j) .eq. 1.0) then 
a (i, j) = a (i, j) + ss 
end if 
end do 
end do 

enddo 


c now find min and max and rescale 



c call f fmax (amax, amin, imax, jmax, imin, jmin, a, isz) 

c amm= ama x - ami n 

do j=l,isz 
do i=l,isz 

c a (i, j ) = (a (i, j ) -amin) /amm 

if (a (i, j) . It . 0 . ) a (i, j) =0 . 0 
enddo 
enddo 

return 

end 

c 

Q ************************************************************* 

c 

function rmsdiff (al, a2, isz) 
real al(isz,isz), a2(isz,isz) 
f sum = 0.0 
do j = 1, isz 

do i = 1, isz 

fsum = fsum + (al(i,j) - a2(i,j))**2 
al (i, j) = a2 (i, j) 
end do 
end do 

rmsdiff = sqrt (fsum) 
return 
end 
c 

0 ***************************************************** 

C 

subroutine rolloff (a, cf, isz) 
complex a (isz, isz) 
real cf 
iszl=isz/2 
do j=l,iszl 
do i=l,iszl 

dl=sqrt ( (i-1. ) **2+( j-1.) **2) /cf 
d2=sqrt ( (i-1. ) **2+ ( j) **2) /cf 
d3=sqrt ( (i) **2+ ( j-1 . ) **2) /cf 
d4=sqrt ( (l.*i) **2+(l . *j) **2) /cf 
a(i, j ) =a (i, j) *cwind (dl) 
a(i,isz+l-j)=a(i,isz+l-j) *cwind (d2 ) 
a (isz+l-i, j) =a (isz+l-i, j) *cwind(d3) 
a (isz+l-i, isz+1- j) =a (isz+l-i, is z+l-j) *cwind(d4) 
enddo 
enddo 
return 
end 
c 

0 ******************************************************************** 
c 
c 
c 

subroutine scale (a, power, isz) 
real a (isz, isz), power 
real sum 
sum=0 . 0 
do j=l,isz 
do i=l,isz 

sum = sum + abs(a(i,j)) 
enddo 
enddo 

rsc = power/sum 
do j=l, isz 
do i=l,isz 



a (i, j) = a (i, j) *rsc 
enddo 
enddo 
return 
end 


c 

Q ****★★*★***★★★★*★**★**★********★★*****★★★*★**★★★**★***★*★*★★*★***★ 

C 

subroutine shift (ph, nar) 
c 

c shifts an image from origin to center by negating every other 

c frequency 

complex ph (nar, nar) 
c 

if (mod(nar,2) .eq.0)then 
nl=nar-l 
n2=nar 
else 
nl=nar 
n2=nar-l 
endif 

do j=2,n2,2 
do i=l,nl,2 

ph ( i, j ) =-ph (i, j ) 
ph ( j, i) =-ph ( j, i) 

enddo 

enddo 

c 

return 

end 

c 

Q ★**★*****★★***★★*★**★★★**************★**★★*★★**★★***★ 

subroutine sqmask (ob j , iradm, isz) 
real obj(isz,isz) 
integer iradm 


do j=l,isz 

do i=l, iradm 

clc= (1 . *iradm-i) / (1 . *iradm) 
obj (i, j) =obj (i, j) *cwind(clc) 
obj ( j, i) =obj ( j, i) *cwind(clc) 
obj (isz-i+1, j) =ob j (isz-i+1, j ) *cwind (clc) 
obj ( j , isz-i+1) =ob j ( j, isz-i+1) *cwind (clc) 
enddo 
enddo 
return 
end 

c 

subroutine subtract (ac, acf, ksz) 
complex ac (ksz, ksz) , acf (ksz, ksz) 
do j=l,ksz 
do i=l,ksz 

ac (i, j) =ac (i, j) -acf (i, j) 
enddo 
enddo 

c ac (1, 1) =cabs (acf (1, 1) ) 

return 



end 


c 


subroutine support <a, mask, test, testn, ksz2) 
c Apply support only to psf 

real a(ksz2), mask(ksz2) 

sump=0 . 

sumn=0 . 

npos=0 

nneg=0 

do j=l,ksz2 

a ( j)=a ( j) *mask ( j) 
if (a(j) .gt.0.)then 
sump=sump+a ( j ) 
npos=npos+l 
endif 

if (a(j) .lt.0.)then 
sumn=sumn+a ( j ) 
nneg=nneg+l 
endif 
enddo 

test=abs (sumn/sump) 
testn=float (nneg) /npos 

return 

end 


c 


c 

c 


★★★★★★★★★A****************************************** 

subroutine wienerd (sig, ref, isz, cf , wl) 
complex sig (isz, isz) 
complex ref (isz, isz) 

iszl=isz/2 

1st calc const level for wiener minimum 
icf=cf 

if(icf.gt.iszl) icf=iszl 
at=0 . 0 
nn=0 

do 100 j=l,iszl 
do 100 i=l,iszl 

idl=sqrt ( (i-1 . ) **2 + ( j-1 . ) **2 ) +0 . 5 
id2=sqrt ( (i-1. ) **2+( j) **2) +0.5 
id3=sqrt ( (i) **2+ < j-1 . ) **2) +0 . 5 
id4=sqrt ( (1 ,*i) **2+(l . *j) **2) +0.5 
if (idl .eq. icf ) then 

at = at+cabs (ref (i, j) ) 
nn=nn+l 
endif 

if (id2 .eq. icf ) then 

at = at+cabs (ref (i, isz+l-j) ) 
nn=nn+l 
endif 

if (id3 . eq. icf ) then 

at = at+cabs (ref (isz+l-i, j) ) 
nn=nn+l 
endif 

if (id4 . eq. icf ) then 

at = at+cabs (ref (isz+l-i, isz+l-j) ) 
nn=nn+l 
endif 
continue 


100 



c 

at=wl*at/nn 

print *, ' Wiener filter constant ' , at 
at2=at*at 
c 

c perform wiener division 
do 200 j=l,isz 
do 200 i=l,isz 

bbl=cabs (ref (i, j) ) **2+at2 
if (bbl . ne . 0 . 0 ) then 

sig (i, j) =sig (i, j) *conjg(ref (i, j) ) / (bbl) 
else 

sig (i, j ) =0.0 
end if 

200 continue 

c 

return 

end 

c 

subroutine wiener f (c, g, f , isz, cf , wl, beta) 
c wiener filter for BID. Adds fraction of previous iteration to 
c wiener filtered result 

complex c (isz, isz) 

complex g (isz, isz) , f (isz, isz) 

real beta,wl,cf 

iszl=isz/2 

c 1st calc const level for wiener minimum 
icf=cf 

if (icf .gt . iszl) icf=iszl 

at=0 . 0 

nn=0 

do j=l,iszl 
do i=l,iszl 

idl=sqrt ( (i-1. ) **2+ ( j-1. ) **2) +0.5 
id2=sqrt ( (i-1. ) **2+( j) **2) +0.5 
id3=sqrt ( (i) **2+( j-1.) **2) +0.5 
id4=sqrt ( (1 . *i) **2+ (1 . * j) **2) +0 . 5 
if (idl .eq. icf ) then 
at = at+cabs (g (i, j ) ) 
nn=nn+l 
endif 

if (id2 .eq. icf ) then 

at = at+cabs (g (i, isz+l-j) ) 
nn=nn+l 
endif 

if (id3 .eq. icf ) then 

at = at+cabs (g (isz+l-i, j) ) 
nn=nn+l 
endif 

if (id4 .eq. icf ) then 

at = at+cabs (g (isz+l-i, isz+l-j) ) 
nn=nn+l 
endif 
enddo 
enddo 


at=wl*at/nn 

at2=at*at 



c perform wiener division 
do j=l,isz 
do i=l,isz 

bbl=cabs (g (i, j) ) **2+at2 
if (bbl . ne . 0 . 0 ) then 

f (i/ j ) = (1 • 0-beta) *f (i, j) +beta*c (i, j) *con jg (g (i, j ) ) / (bbl) 
end if 
enddo 
enddo 

return 

end 

c 

Q************************************ ********************** 

C 

c 

subroutine zero(a,isz) 
complex a(isz,isz) 
do j=l,isz 
do i=l,isz 
a (i, j) =0 . 0 
enddo 
enddo 
return 
end 
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Fourt . f 

NN, ISIGN AND IFORM MUST ALL 
BE DIMENSIONED INTEGER* 4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

THE FAST FOURIER TRANSFORM IN USAS I BASIC FORTRAN 

TRANSFORM ( Jl, J2, . . . ) = SUM(DATA(I1, 12, . . . ) *W1** ( (11-1) * (Jl-1) ) 

*W2** ( (12-1) * (J2-1) )*...) , 

WHERE II AND Jl RUN FROM 1 TO NN(1) AND W1=EXP ( ISIGN*2*PI * 

SQRT ( - 1 ) / NN ( 1 ) ) , ETC. THERE IS NO LIMIT ON THE DIMENSIONALITY 
(NUMBER OF SUBSCRIPTS) OF THE DATA ARRAY. IF AN INVERSE 
TRANSFORM (ISIGN=+1) IS PERFORMED UPON AN ARRAY OF TRANSFORMED 
(ISIGN=-1) DATA, THE ORIGINAL DATA WILL REAPPEAR, 

MULTIPLIED BY NN ( 1 ) *NN (2 ) * . . . THE ARRAY OF INPUT DATA MAY BE 
REAL OR COMPLEX, AT THE PROGRAMMERS OPTION, WITH A SAVING OF 
UP TO FORTY PER CENT IN RUNNING TIME FOR REAL OVER COMPLEX. 

(FOR FASTEST TRANSFORM OF REAL DATA, NN(1) SHOULD BE EVEN.) 

THE TRANSFORM VALUES ARE ALWAYS COMPLEX, AND ARE RETURNED IN THE 
ORIGINAL ARRAY OF DATA, REPLACING THE INPUT DATA. THE LENGTH 
OF EACH DIMENSION OF THE DATA ARRAY MAY BE ANY INTEGER. THE 
PROGRAM RUNS FASTER ON COMPOSITE INTEGERS THAN ON PRIMES, AND IS 
PARTICULARLY FAST ON NUMBERS RICH IN FACTORS OF TWO. 

TIMING IS IN FACT GIVEN BY THE FOLLOWING FORMULA. LET NTOT BE THE 
TOTAL NUMBER OF POINTS (REAL OR COMPLEX) IN THE DATA ARRAY, THAT 
IS, NTOT=NN ( 1 ) *NN ( 2 ) * . . . DECOMPOSE NTOT INTO ITS PRIME FACTORS, 
SUCH AS 2**K2 * 3**K3 * 5**K5 * ... LET SUM2 BE THE SUM OF ALL 
THE FACTORS OF TWO IN NTOT, THAT IS, SUM2 = 2*K2. LET SUMF BE 
THE SUM OF ALL OTHER FACTORS OF NTOT, THAT IS, SUMF = 3*K3+5*K5+.. 
THE TIME TAKEN BY A MULTIDIMENSIONAL TRANSFORM ON THESE NTOT DATA 
POINT ADD TIME = SIX MICROSECONDS), T = 3000 + NTOT* ( 600+40*SUM2+ 
IS T = TO + NTOT* (T1+T2*SUM2+T3*SUMF) . ON THE CDC 3300 (FLOATING 
175*SUMF) MICROSECONDS ON COMPLEX DATA. 

IMPLEMENTATION OF THE DEFINITION BY SUMMATION WILL RUN IN A TIME 
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PROPORTIONAL TO NTOT* (NN ( 1 ) +NN (2 )+...) . FOR HIGHLY COMPOSITE NTOT 
THE SAVINGS OFFERED BY THIS PROGRAM CAN BE DRAMATIC. A ONE-DIMEN- 
SIONAL ARRAY 4000 IN LENGTH WILL BE TRANSFORMED IN 4000* (600+ 

40* (2+2+2+2+2) +175* (5+5+5) ) = 14.5 SECONDS VERSUS ABOUT 4000* 
4000*175 = 2800 SECONDS FOR THE STRAIGHTFORWARD TECHNIQUE. 

THE FAST FOURIER ALGORITHM PLACES TWO RESTRICTIONS UPON THE 

NATURE OF THE DATA BEYOND THE USUAL RESTRICTION THAT 

THE DATA FORM ONE CYCLE OF A PERIODIC FUNCTION. THEY ARE — 

1. r THE NUMBER OF INPUT DATA AND THE NUMBER OF TRANSFORM VALUES 
MUST BE THE SAME. 

2. CONSIDERING THE DATA TO BE IN THE TIME DOMAIN, 

THEY MUST BE EQUI-SPACED AT INTERVALS OF DT. FURTHER, THE TRANS- 
FORM VALUES, CONSIDERED TO BE IN FREQUENCY SPACE, WILL BE EQUI- 
SPACED FROM 0 TO 2 *PI* (NN ( I ) -1 ) / (NN ( I ) *DT) AT INTERVALS OF 
2*PI/ (NN (I) *DT) FOR EACH DIMENSION OF LENGTH NN(I). OF COURSE, 

DT NEED NOT BE THE SAME FOR EVERY DIMENSION. 

THE CALLING SEQUENCE IS — 

CALL FOURT (DATA, NN, NDIM, ISIGN, IFORM, WORK) 

DATA IS THE ARRAY USED TO HOLD THE REAL AND IMAGINARY PARTS 
OF THE DATA ON INPUT AND THE TRANSFORM VALUES ON OUTPUT. IT 
IS A MULTIDIMENSIONAL FLOATING POINT ARRAY, WITH THE REAL AND 
IMAGINARY PARTS OF A DATUM STORED IMMEDIATELY ADJACENT IN STORAGE 
(SUCH AS FORTRAN IV PLACES THEM) . THE EXTENT OF EACH DIMENSION 
IS GIVEN IN THE INTEGER ARRAY NN, OF LENGTH NDIM. ISIGN IS -1 
TO INDICATE A FORWARD TRANSFORM (EXPONENTIAL SIGN IS -) AND +1 
FOR AN INVERSE TRANSFORM (SIGN IS +) . IFORM IS +1 IF THE DATA AND 
THE TRANSFORM VALUES ARE COMPLEX. IT IS 0 IF THE DATA ARE REAL 
BUT THE TRANSFORM VALUES ARE COMPLEX. IF IT IS 0, THE IMAGINARY 
PARTS OF THE DATA SHOULD BE SET TO ZERO. AS EXPLAINED ABOVE, THE 
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TRANSFORM VALUES ARE ALWAYS COMPLEX AND ARE STORED IN ARRAY DATA. 
WORK IS AN ARRAY USED FOR WORKING STORAGE. IT IS NOT NECESSARY 
IF ALL THE DIMENSIONS OF THE DATA ARE POWERS OF TWO. IN THIS CASE 
IT MAY BE REPLACED BY 0 IN THE CALLING SEQUENCE. THUS, USE OF 
POWERS OF TWO CAN FREE A GOOD DEAL OF STORAGE. IF ANY DIMENSION 
IS NOT A POWER OF TWO, THIS ARRAY MUST BE SUPPLIED. IT IS 
FLOATING POINT, ONE DIMENSIONAL OF LENGTH EQUAL TO TWICE THE 
LARGEST ARRAY DIMENSION (I.E., NN(I) ) THAT IS NOT A POWER OF 
TWO. THEREFORE, IN ONE DIMENSION FOR A NON POWER OF TWO, 

WORK OCCUPIES AS MANY STORAGE LOCATIONS AS DATA. IF SUPPLIED, 

WORK MUST NOT BE THE SAME ARRAY AS DATA. ALL SUBSCRIPTS OF ALL 
ARRAYS BEGIN AT ONE. 

EXAMPLE 1. THREE-DIMENSIONAL FORWARD FOURIER TRANSFORM OF A 
COMPLEX ARRAY DIMENSIONED 32 BY 25 BY 13 IN FORTRAN IV. 

COMPLEX ARRAY DIMENSIONED 32 BY 25 BY 13 IN FORTRAN IV. 

DIMENSION DATA (32, 25, 13) ,WORK(50) ,NN(3) 

COMPLEX DATA 
DATA NN/32, 25, 13/ 

DO 1 1=1,32 
DO 1 J=l, 25 
DO 1 K=l, 13 

DATA ( I , J, K) =COMPLEX VALUE 
CALL FOURT (DATA,NN, 3, -1, 1, WORK) 

EXAMPLE 2. ONE -DIMENSIONAL FORWARD TRANSFORM OF A REAL ARRAY OF 
LENGTH 64 IN FORTRAN II. 

DIMENSION DATA ( 2 , 64 ) 

DO 2 1=1, 64 

DATA ( 1 , I ) =REAL P ART 

DATA (2, I) =0 . 

CALL FOURT (DATA, 64, 1,-1, 0,0) 
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# 126 
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THERE ARE NO ERROR MESSAGES OR ERROR HALTS IN THIS PROGRAM. THE 
PROGRAM RETURNS IMMEDIATELY IF NDIM OR ANY NN(I) IS LESS THAN ONE. 

PROGRAM BY NORMAN BRENNER FROM THE BASIC PROGRAM BY CHARLES 
RADER (BOTH OF MIT LINCOLN LABORATORY). MAY 1967. THE IDEA 
FOR THE DIGIT REVERSAL WAS SUGGESTED BY RALPH ALTER (ALSO MIT LL) . 
THIS IS THE FASTEST AND MOST VERSATILE VERSION OF THE FFT KNOWN 
TO THE AUTHOR. A PROGRAM CALLED FOUR2 IS AVAILABLE THAT ALSO 
PERFORMS THE FAST FOURIER TRANSFORM AND IS WRITTEN IN USASI BASIC 
FORTRAN. IT IS ABOUT ONE THIRD AS LONG AND RESTRICTS THE 
DIMENSIONS OF THE INPUT ARRAY (WHICH MUST BE COMPLEX) TO BE POWERS 
OF TWO. ANOTHER PROGRAM, CALLED FOUR1, IS ONE TENTH AS LONG AND 
RUNS TWO THIRDS AS FAST ON A ONE -DIMENSIONAL COMPLEX ARRAY WHOSE 
LENGTH IS A POWER OF TWO. 

REFERENCE — 

FAST FOURIER TRANSFORMS FOR FUN AND PROFIT, W. GENTLEMAN AND 
G. SANDE, 1966 FALL JOINT COMPUTER CONFERENCE. 

THE WORK REPORTED IN THIS DOCUMENT WAS PERFORMED AT LINCOLN LAB- 
ORATORY, A CENTER FOR RESEARCH OPERATED BY MASSACHUSETTS INSTITUTE 
OF TECHNOLOGY, WITH THE SUPPORT OF THE U.S. AIR FORCE UNDER 
CONTRACT AF 1 9 ( 628 ) -51 67 . 

subroutine fourt(data, nn, ndim, isign, iform, work) 
dimension data(l), nn(l), if act (32), work(l) 

"fourt . for" 
twopi = 6.283185307 
rthlf = .7071067812 
if (ndim - 1) 920, 1, 1 
ntot = 2 

do 2 idim = 1, ndim 
if (nn(idim)) 920, 920, 2 


MAIN LOOP FOR EACH DIMENSION 
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"fourt . for" 

ntot — ntot * nn(idim) 

"fourt . for" 

npl = 2 

do 910 idim = 1, ndim 
n = nn(idim) 
np2 = npl * n 


IS N A POWER OF TWO AND IF NOT, WHAT ARE ITS FACTORS 


"fourt. for" 

if (n - 1) 920, 900, 5 

"fourt . for" 

m = n 

ntwo = npl 

if = 1 

idiv = 2 

iquot = m / idiv 

irem = m - (idiv * iquot) 

if (iquot - idiv) 50, 11, 11 

if (irem) 20, 12, 20 

ntwo = ntwo + ntwo 

ifact(if) = idiv 

if = if + 1 

m = iquot 

goto 10 

idiv = 3 

inon2 = if 

iquot = m / idiv 

irem = m - (idiv * iquot) 

if (iquot - idiv) 60, 31, 31 

if (irem) 40, 32, 40 

ifact(if) = idiv 

if = if + 1 

m = iquot 

goto 30 

idiv = idiv + 2 
goto 30 
inon2 = if 

if (irem) 60, 51, 60 
ntwo = ntwo + ntwo 
goto 70 
ifact (if) = m 


SEPARATE FOUR CASES — 

1. COMPLEX TRANSFORM 

2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION. METHOD 
TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CON 
JUGATE SYMMETRY. 

3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD — 
SET THE IMAGINARY PARTS TO ZERO. 

4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD — 



TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS 


c 
c 

C ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS 

C 

C ARE THE ODD NUMBERED REAL VALUES . SEPARATE AND SUPPLY 

C 

C THE SECOND HALF BY CONJUGATE SYMMETRY. 

C 

C 

C 

# 174 "f ourt . for" 

70 non2p = np2 / ntwo 

# 189 "fourt . for" 

lease = 1 
ifmin = 1 
ilrng = npl 

if (idim - 4) 74, 100, 100 
74 if (iform) 71, 71, 100 

71 icase = 2 

ilrng = npO * (1 + (nprev / 2) ) 
if (idim - 1) 72, 72, 100 

72 icase = 3 
ilrng = npl 

if (ntwo - npl) 100, 100, 73 

73 icase = 4 
ifmin = 2 

ntwo = ntwo / 2 
n = n / 2 
np2 = np2 / 2 
ntot = ntot / 2 
i = 1 

do 80 j = 1, ntot 
data(j) = data(i) 
c 
c 

C SHUFFLE DATA BY BIT REVERSAL, SINCE N=2**K. AS THE SHUFFLING 

C 

C CAN BE DONE BY SIMPLE INTERCHANGE, NO WORKING ARRAY IS NEEDED 

C 
C 
C 

# 209 "fourt. for" 

80 i = i + 2 

# 214 "fourt. for" 

100 if (non2p - 1) 101, 101, 200 

101 np2hf = np2 / 2 

j = 1 

do 150 i2 = 1, np2, npl 
if (j - i2 ) 121, 130, 130 
121 ilmax = (i2 + npl) - 2 
do 125 il = i2, ilmax, 2 
do 125 i3 = il, ntot, np2 
j3 = (j + i3) - i2 
tempr = data(i3) 
tempi = data(i3 + 1) 
data(i3) = data(j3) 
data(i3 +1) = data(j3 + 1) 
data(j3) = tempr 
125 data(j3 +1) = tempi 
130 m = np2hf 

140 if (j - m) 150, 150, 141 

141 j = j - m 
m = m / 2 

if (m - npl) 150, 140, 140 



150 j = j + m 
c 
c 

C SHUFFLE DATA BY DIGIT REVERSAL FOR GENERAL N 

C 
C 
C 

# 235 "f ourt . for" 

goto 300 

# 239 "f ourt. for" 

200 nwork = 2 * n 

do 270 il = 1, npl, 2 
do 270 13 = il, ntot, np2 
j = 13 

do 260 1=1, nwork, 2 
if (icase - 3) 210, 220, 210 
210 work (i) = data(j) 

work(i +1) = data(j + 1) 
goto 240 

220 work(i) = data(j) 
work (i +1) =0. 

240 ifp2 = np2 
if = ifmin 

250 ifpl = ifp2 / ifact(if) 
j = j + ifpl 

if ((j - i3) - ifp2 ) 260, 255, 255 
255 j = j - ifp2 
ifp2 = ifpl 
if = if + 1 

if (ifp2 - npl) 260, 260, 250 
260 continue 

i2max = (i3 + np2) - npl 
i = 1 

do 270 i2 = i3, i2max, npl 
data(i2) = work(i) 
data(i2 + 1) = work(i + 1) 

C MAIN LOOP FOR FACTORS OF TWO. 

C 

c W=EXP (ISIGN*2*PI*SQRT (-1) *M/ (4*MMAX) ) . CHECK FOR W=ISIGN*SQRT ( -1 ) 

C 

C AND REPEAT FOR W=W* ( 1+ISIGN*SQRT (-1) ) /SQRT (2) . 

C 

C 

C 

# 265 "fourt.for" 

270 i = i + 2 

# 270 "fourt.for" 

300 if (ntwo - npl) 600, 600, 305 
305 npltw = npl + npl 
ipar = ntwo / npl 
310 if (ipar - 2) 350, 330, 320 
320 ipar = ipar / 4 
goto 310 

330 do 340 il = 1, ilrng, 2 

do 340 kl = il, ntot, npltw 
k2 = kl + npl 
tempr = data(k2) 
tempi = data(k2 + 1) 
data(k2) = data(kl) - tempr 
data (k2 + 1) = data (kl + 1) - tempi 
data(kl) = data(kl) + tempr 
340 data(kl + 1) = data (kl + 1) + tempi 
350 mmax = npl 

360 if (mmax - (ntwo / 2) ) 370, 600, 600 


370 lmax = maxO (npltw,mmax / 2) 
do 570 1 = npl, lmax, npltw 
m = 1 

if (mmax - npl) 420, 420, 380 
380 theta = - ( (twopi * float (1)) / float (4 * mmax)) 
if (isign) 400, 390, 390 
390 theta = - theta 
400 wr = cos (theta) 
wi = sin (theta) 

410 w2r = (wr * wr) - (wi * wi) 
w2i = (2. * wr) * wi 
w3r = (w2r * wr) - (w2i * wi) 
w3i = (w2r * wi) + (w2i * wr) 

420 do 530 il = 1, ilrng, 2 
kmin = il + (ipar * m) 
if (mmax - npl) 430, 430, 440 
430 kmin = il 
440 kdif = ipar * mmax 
450 kstep = 4 * kdif 

if (kstep - ntwo) 460, 460, 530 
460 do 520 kl = kmin, ntot, kstep 

k2 = kl + kdif 

k3 = k2 + kdif 

k4 = k3 + kdif 

if (mmax - npl) 470, 470, 480 

470 ulr = data(kl) + data(k2) 

uli = data (kl + 1) + data(k2 + 1) 

u2r = data(k3) + data(k4) 

u2i = data(k3 + 1) + data(k4 + 1) 

u3r = data(kl) - data(k2) 

u3i = data (kl + 1) - data(k2 + 1) 

if (isign) 471, 472, 472 

471 u4r = data(k3 + 1) - data (k4 + 1) 
u4i = data(k4) - data(k3) 

goto 510 

472 u4r = data(k4 + 1) - data(k3 + 1) 
u4i = data(k3) - data(k4) 

goto 510 

480 t2r = (w2r * data(k2)) - (w2i * data(k2 + 1)) 
t2i = (w2r * data(k2 +1)) + (w2i * data(k2)) 
t3r = (wr * data(k3)) - (wi * data (k3 + 1)) 

t3i = (wr * data(k3 + 1) ) + (wi * data (k3) ) 

t4r = (w3r * data(k4)) - (w3i * data (k4 + 1)) 
t4i = (w3r * data(k4 + 1) ) + (w3i * data(k4)) 
ulr = data(kl) + t2r 
uli = data(kl + 1) + t2i 
u2r = t3r + t4r 

u2i = t3i + t4i 

u3r = data(kl) - t2r 
u3i = data(kl + 1) - t2i 
if (isign) 490, 500, 500 
490 u4r = t3i - t4i 

u4i = t4r - t3r 

goto 510 

500 u4r = t4i - t3i 

u4i = t3r - t4r 

510 data(kl) = ulr + u2r 

data(kl + 1) = uli + u2i 
data(k2) = u3r + u4r 
data(k2 + 1) = u3i + u4i 
data(k3) = ulr - u2r 
data(k3 + 1) = uli - u2i 
data(k4) = u3r - u4r 
520 data(k4 + 1) = u3i - u4i 



kdif = kstep 

kmin = (4 * (kmin - il) ) + il 
goto 450 
530 continue 

m = m + lmax 

if (m - max) 540, 540, 570 
540 if (isign) 550, 560, 560 
550 tempr = wr 

wr = (wr + wi) * rthlf 
wi = (wi - tempr) * rthlf 
goto 410 
560 tempr = wr 

wr = (wr - wi) * rthlf 
wi = (tempr + wi) * rthlf 
goto 410 
570 continue 

ipar = 3 - ipar 
max = max + mmax 
c 
c 

C MAIN LOOP FOR FACTORS NOT EQUAL TO TWO. 

C 

C W=EXP (ISIGN*2*PI*SQRT (-1) * (J1+J2-I3-1) /IFP2) 

C 

c 

c 

# 369 "fourt. for" 

goto 360 

# 374 "fourt. for" 

600 if (non2p - 1) 700, 700, 601 

601 ifpl = ntwo 
if = inon2 

610 ifp2 = ifact(if) * ifpl 

theta - - (twopi / float (ifact (if )) ) 
if (isign) 612, 611, 611 

611 theta = - theta 

612 wstpr = cos (theta) 
wstpi = sin (theta) 

do 650 jl = 1, ifpl, npl 

thetm = - ((twopi * float (jl -1)) / float (ifp2) ) 
if (isign) 614, 613, 613 

613 thetm = - thetm 

614 wminr = cos (thetm) 
wmini = sin (thetm) 
ilmax = (jl + ilrng) - 2 
do 650 il = jl, ilmax, 2 
do 650 i3 = il, ntot, np2 
i = 1 

wr = wminr 
wi = wmini 

j2max = (i3 + ifp2) - ifpl 
do 640 j2 = i3, j2max, ifpl 
twowr = wr + wr 
j3max = (j2 + np2) - ifp2 
do 630 j3 = j2, j3max, ifp2 
jmin = (j3 - j2) + i3 
j = (jmin + ifp2) - ifpl 
sr = data(j) 
si = data(j + 1) 
oldsr = 0. 
oldsi = 0. 
j = j - ifpl 
620 stmpr = sr 
stmpi = si 



sr = ( (twowr * sr) - oldsr) + data(j) 
si = ({twowr * si) - oldsi) + data(j + 1) 
oldsr = stmpr 
oldsi = stmpi 
j = j - ifpl 

if (j - jmin) 621, 621, 620 

621 work(i) = ( ( (wr * sr) - (wi * si)) - oldsr) + data(j) 

work(i + 1) = ( ( (wi * sr) + (wr * si)) - oldsi) + data(j + 1) 
630 i = i + 2 

wtemp = wr * wstpi 

wr = (wr * wstpr) - (wi * wstpi) 

640 wi = (wi * wstpr) + wtemp 
i = 1 

do 650 j2 = i3, j2max, ifpl 
j3max = (j2 + np2) - ifp2 
do 650 j3 = j2, j3max, ifp2 
data(j3) = work(i) 
data(j3 +1) = work(i + 1) 

650 i = i + 2 
if = if + 1 
ifpl = ifp2 
c 
c 

C COMPLETE A REAL TRANSFORM IN THE 1ST DIMENSION, N EVEN, BY CON- 

C 

C JUGATE SYMMETRIES . 

c 
c 
c 

# 430 "fourt.for" 

if (ifpl - np2) 610, 700, 700 

# 435 "fourt.for" 

700 goto (900, 800, 900, 701), icase 

701 nhalf = n 
n = n + n 

theta = - (twopi / float (n) ) 
if (isign) 703, 702, 702 

702 theta = - theta 

703 wstpr = cos (theta) 
wstpi = sin (theta) 
wr = wstpr 

wi = wstpi 
imin = 3 

jmin = (2 * nhalf) ' - 1 
goto 725 
710 j = jmin 

do 720 i = imin, ntot, np2 
sumr = (data(i) + data(j)) / 2. 
sumi = (data(i + 1) + data(j + 1)) / 2. 
difr = (data(i) - data(j)) / 2. 
difi = (data(i + 1) - data(j +1)) / 2. 
tempr = (wr * sumi) + (wi * difr) 

tempi = (wi * sumi) - (wr * difr) 

data(i) = sumr + tempr 
data(i + 1) = difi + tempi 
data(j) = sumr - tempr 
data(j + 1) = (- difi) + tempi 
720 j = j + np2 

imin = imin + 2 

jmin = jmin - 2 

wtemp = wr * wstpi 
wr = (wr * wstpr) - (wi * wstpi) 
wi = (wi * wstpr) + wtemp 
725 if (imin - jmin) 710, 730, 740 



730 if (isign) 731, 740, 740 

731 do 735 i = imin, ntot, np2 
735 data(i +1) = - data(i + 1) 

740 np2 = np2 + np2 

ntot = ntot + ntot 
j = ntot + 1 
imax = (ntot / 2) + 1 
745 imin = imax - (2 * nhalf) 
i = imin 
goto 755 

750 data(j) = data(i) 

data(j + 1) = - data(i + 1) 

755 i = i + 2 

j = j - 2 

if (i - imax) 750, 760, 760 
760 data(j) = data (imin) - data (imin + 1) 
data(j +1) =0. 
if (i - j) 770, 780, 780 
765 data(j) = data(i) 

data(j + 1) = data(i + 1) 

770 i = i - 2 

j = j - 2 

if (i - imin) 775, 775, 765 
775 data(j) = data (imin) + data (imin + 1) 
data(j +1) =0. 
imax = imin 
goto 745 

780 data (1) = data(l) + data (2) 
data ( 2 ) = 0 . 
c 
c 

C COMPLETE A REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION BY 

C 

C CONJUGATE SYMMETRIES . 

c 
c 
c 

# 496 "fourt.for" 

goto 900 

# 501 "fourt.for" 

800 if (ilrng - npl) 805, 900, 900 
805 do 860 i3 = 1, ntot, np2 
i2max = (i3 + np2) - npl 
do 860 i2 = i3, i2max, npl 
imax = (i2 + npl) - 2 
imin = i2 + ilrng 
jmax = ( (2 * i3) + npl) - imin 
if (i2 - i3) 820, 820, 810 
810 jmax = jmax + np2 
820 if (idim - 2) 850, 850, 830 
830 j = jmax + npO 

do 840 i = imin, imax, 2 

data (i) = data ( j) 

data(i + 1) = - data(j + 1) 

840 j = j - 2 
850 j = jmax 

do 860 i = imin, imax, npO 
data(i) = data(j) 
data(i +1) = - data(j + 1) 
c 
c 

C END OF LOOP ON EACH DIMENSION 

C 

C 



c 

# 520 "fourt.for" 
860 j = j - npO 

# 524 "fourt.for" 
900 npO = npl 

npl = np2 
910 nprev = n 
920 return 
end 



